rm(list=ls())

##############################################################
#### 0. Load Packages 
##############################################################

packages <- c("ggplot2", "cshapes", "haven", "dplyr", "RColorBrewer")

package_installed <-
  sapply(packages, function(pack)
    pack %in% rownames(installed.packages()))

if (any(!package_installed)) {
  sapply(packages[!package_installed], install.packages)
}

sapply(packages, require, character.only = TRUE)

rm(packages,package_installed)

##############################################################
#### 1. Load Data 
##############################################################

trac <- read_dta("COVID19_TRAC.dta")

names(trac)[names(trac) == "cowcode"] <- "COWCODE"
trac$COWCODE <- ifelse(trac$COWCODE == "1015", "345", trac$COWCODE)
trac$COWCODE <- ifelse(trac$COWCODE == "817", "816", trac$COWCODE)
trac$COWCODE <- as.integer(trac$COWCODE)

trac$be_china <- ifelse(trac$be_cat_1_China_carryover==1, 3, ifelse(trac$be_cat_2_China_carryover==1, 2, ifelse(trac$be_cat_3_China_carryover==1, 1, 0)))
trac$be_china <- factor(trac$be_china, levels = c(0,1,2,3), labels = c("No restrictions", "Category 3", "Category 2", "Category 1"))

##############################################################
#### 2-a. Figure 2-a: Entry Restrictions against China (2020-02-15)
##############################################################

trac_subset <- subset(trac, date == "2020-02-15")

world <- cshp(date=as.Date("2016-1-1"))
world.df <- fortify(world)
world.df <- data.frame(world.df, world@data[world.df$id, ])
world.df.subset <- left_join(world.df, trac_subset, by=c("COWCODE"), all = TRUE)

pdf("Figure2a.pdf", width=9, height=5)
mapPlot <- ggplot() + theme_bw() + theme(plot.background = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.border = element_blank(),axis.title = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),legend.position="bottom",legend.key.width = unit(1, "cm"))
mapPlot <- mapPlot + geom_polygon(data=world.df.subset, aes(x = long, y = lat, group = group, fill = be_china), lwd = 1/5, colour="gray60") + scale_fill_manual(name = "", values = c("white", "#FDBB84", "#EF6548", "#7F0000"), drop=FALSE, na.translate = FALSE) + theme(legend.key.width = unit(1, "cm")) + guides(col = guide_legend(nrow = 1, reverse = TRUE)) + theme(legend.title=element_text(family="Times", size=17), legend.text = element_text(family="Times", size = 20)) 
plot(mapPlot)
dev.off()

##############################################################
#### 2-b. Figure 2-b: Entry Restrictions against China (2020-03-15)
##############################################################

trac_subset <- subset(trac, date == "2020-03-15")

world <- cshp(date=as.Date("2016-1-1"))
world.df <- fortify(world)
world.df <- data.frame(world.df, world@data[world.df$id, ])
world.df.subset <- left_join(world.df, trac_subset, by=c("COWCODE"), all = TRUE)

pdf("Figure2b.pdf", width=9, height=5)
mapPlot <- ggplot() + theme_bw() + theme(plot.background = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.border = element_blank(),axis.title = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),legend.position="bottom",legend.key.width = unit(1, "cm"))
mapPlot <- mapPlot + geom_polygon(data=world.df.subset, aes(x = long, y = lat, group = group, fill = be_china), lwd = 1/5, colour="gray60") + scale_fill_manual(name = "", values = c("white", "#FDBB84", "#EF6548", "#7F0000"), drop=FALSE, na.translate = FALSE) + theme(legend.key.width = unit(1, "cm")) + guides(col = guide_legend(nrow = 1, reverse = TRUE)) + theme(legend.title=element_text(family="Times", size=17), legend.text = element_text(family="Times", size = 20)) 
plot(mapPlot)
dev.off()

##############################################################
#### 4. Figure 4: TRAC and Passenger Volume
##############################################################

load("oag_flight.RData")

ggplot(subset(oag, month > 201906 & month < 202010), aes(x = calendar_month, y = passenger_china, fill=as.factor(category))) + 
  geom_bar(stat="identity") + theme_minimal() + facet_wrap(~country, nrow = 4) + theme(legend.position="bottom") + theme(legend.title = element_blank()) +
  scale_fill_manual(values=c("#999999", "#F8766D", "#00BA38", "#619CFF")) +
  ylab("Flight Passengers from China") + xlab("Month")

ggsave("Figure4.pdf", width = 16, height = 18, unit = "cm")

